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Adaptation to environmental change is a common property of biological systems. Cells initially 
respond to external changes in the environment, but after some time, they regain their original 
state. By considering an element consisting of two variables that show such adaptation dynamics, 
we studied a coupled dynamical system containing such elements to examine the diverse dynamics 
in the system and classified the behaviors on the basis of the network structure that determined 
the interaction among elements. For a system with two elements, two types of behaviors, perfect 
adaptation and simple oscillation, were observed. For a system with three elements, in addition 
to these two types, novel types of dynamics, namely, rapid burst-type oscillation and a slow cycle, 
were discovered; depending on the initial conditions, these novel types of dynamics coexisted. These 
behaviors are a result of the characteristic dynamics of each element, i.e., fast response and slow 
adaptation processes. The behaviors depend on the network structure (in specific, a combination 
of positive or negative feedback among elements). Cooperativity among elements due to a positive 
feedback loop leads to simple oscillation, whereas frustration involving alternating positive and 
negative interactions among elements leads to the coexistence of rapid bursting oscillation and a 
slow cycle. These behaviors are classified on the basis of the frustration indices defined by the 
network structure. The period of the slow cycle is much longer than the original adaptation time 
scale, while the burst-type oscillation is a continued response that does not involve any adaptation. 
We briefly discuss the universal applicability of our results to a network of a larger number of 
elements and their possible relevance to biological systems. 



PACS numbers: 05.45.-a, 87.10.-e, 82.39. Rt 

I. INTRODUCTION 

Adaptation is a common property of biological organ- 
isms. Individual organisms initially sense and respond 
to external changes in their environment; however, after 
some time, they become habituated to the new environ- 
ment, and hence, the response disappears. Adaptation 
can occur on a time scale of generations, but it also pro- 
gresses on much shorter time scales of seconds, minutes, 
or hours, i.e., within the lifetime of an individual organ- 
ism. Here, we focus on the adaptive behavior on such 
short timescales; such behavior is generally observed not 
only in multicellular organisms but also in unicellular or- 
ganisms like E.coli or Paramecium and is regarded as a 
general property of a cell. 

The study of cellular adaptation behavior has been pi- 
oneered by Koshland[l|, and Oosawa[H H[. Koshland 
classified adaptation as absolute(perfect) and partial 
adaptation. Environmental change has an impact on 
intracellular states such as gene expression patterns 
through reaction networks, for example, by signal trans- 
duction. In absolute adaptation, variables representing 
the concentration of some chemicals start to change in 
response to environmental change, and then they re- 
gain their original values even though the environmental 
change is not reversed. Indeed, several examples of sim- 
ple biochemical circuits have been proposed to represent 



such adaptive behaviors [EH0|- ^ n addition to theoreti- 
cal interest, experimental verification of such adaptation 
processes has gained much attention (§, [oL [ToL IfO] . 

Adaptive response is not limited to a specific variable 
of a biological system. Not one but many variables re- 
gain their original values within a certain time after an 
environmental change has occurred. This type of re- 
sponse has been studied in the microarray analysis of 
gene expression levels [13, EJ- Often, some partial units 
in a whole biological system have a tendency to main- 
tain their state, independent of the external conditions 
— this is known as homeostasis. On the other hand, they 
also have to respond to an external change in order to 
survive. Hence, the adaptive behavior is expected to be 
an inherent process that involves many partial units of 
a biological system. Then, the dynamic behavior of the 
whole system is determined by the behavior of several 
adaptation units that interact with each other. 

Now, in addition to revealing the processes and mech- 
anisms in an adaptation unit, i.e., an adaptation mod- 
ule, it is also important to study a system consisting 
of several adaptation units that are not isolated but in- 
teract with other units. Here, such a unit of adaptive 
behavior can be a cell in a colony of interacting cells 
or in a multiceullar organism. In this case, the interac- 
tion occurs through chemical signals resulting in cell-cell 
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communication. Alternatively, a certain part of a global 
intracellular reaction network can function as a unit to 
exhibit adaptive behavior. In this case, each unit inter- 
acts through reaction paths in the network. For exam- 
ple, a signal transduction network or metabolic network 
is highly complicated, involving many feedback and feed- 
forward interactions or multi step reactions. Each adap- 
tive or non-adaptive unit is embeded within a cell and 
receives inputs not only from outside the cell but also 
from other units within the intracellular network. Such 
interaction between units should play an important role 
in the global function of a cell. 

In this study, we examine the behavior of a coupled 
system consisting of elements that exhibit adaptive be- 
havior. If each clement is a module in a chemical reaction 
network or a reaction motif, the coupled adaptive system 
represents a global intracellular reaction network. If each 
element represents a cell, the system represents the co- 
operative behavior of cells within a multicellular system. 

As for dynamical systems with many degrees of free- 
dom, coupled oscillators 1131 and coupled chaotic systems 
(such as coupled mapsfjlS. [l4| have been extensively 
studied over decades. Several novel concepts such as syn- 
chronization transition, clustering chaotic itinerancy, and 
collective chaos have been developed in these studies. In 
contrast, systems of coupled adaptive elements have not 
been systematically studied thus far. However, such a 
study may help develop novel concepts and reveal non- 
trivial interesting behaviors. In fact, dynamics that facil- 
itate perfect adaptation have the following characteristics 
that distinguish them from other systems: when a spe- 
cific parameter of the system is changed, the response to 
this change appears first as a change in the values of the 
variables, but later, some variables regain their original 
values independent of the parameter values. This inde- 
pendence results in the imposition of some constraints. 
Further, the above behavior implies the existence of two 
time scales: one for response and the other for relaxation 
to the original value. Therefore, the manner of return 
to the original values and the existence of the two time 
scales result in nontrivial dynamics in a coupled system. 
Here, we classify the behaviors of coupled adaptive sys- 
tems on the basis of the interaction network structure 
among elements. For the classification, we will introduce 
the measure of frustration to characterize the interaction 
among elements. 

The paper is organized as follows. In section II, we 
introduce a model of coupled adaptive elements. In sec- 
tion III, we present results for the two-element case. The 
three-element case is studied in section IV. In this case, 
the coexistence of a slow cycle and fast bursting oscilla- 
tions is observed when frustration is taken into account. 
Extensions to cases with a larger number of elements will 
also be discussed. 

II. MODEL 

Here, we introduce a model of coupled adaptive ele- 
ments. In this model, elements that show adaptive be- 



havior to an input interact with each other. Here, adap- 
tation refers to dynamics by which some variables regain 
their original value, independent of the value of the ex- 
ternal inputs. 

Now, we consider a minimal system of such adapta- 
tion. When the value of a constant external input S is 
changed, first, a state variable u changes in response to 
this input, however, after some time, u regains its orig- 
inal value independent of the value of S. This system 
is a type of excitable system in which a state variable 
changes in response to a stimulus and then returns to 
the pre-stimulated state. The simplest way to construct 
such dynamics is to introduce another variable (v) that 
"absorbs" the changes in S. Hence, we consider the fol- 
lowing two- variable dynamical system: 



du 
~dt ' 

dv 

~dt 



f(u,v;S), 
--g(u,v;S), 



(1) 



where the fixed point solution u* , v* given by 
f(u*,v*;S) = 0, g(u*,v*;S) = should be set to be 
stable. The above condition is satisfied if u* is invariant 
against the change in S, even though / increases with S. 
Then, the postulate for the response to S and adapta- 
tion is satisfied because u first increases with S and then 
regains the value u*. 

As an example of such adaptive dynamics, we consider 
a chemical reaction. In the reaction process, time evo- 
lution is governed by a set of rate equations in which 
variables u and v represent chemical concentrations. In- 
deed, there have been several models that exhibit adap- 
tive dynamics @, 0, [IB Gil- In this article, we adopt a 
system obtained by adiabatically eliminating a variable 
from a model originally introduced by Levchenko and 
Iqlesias 17]. In this simplified model, a chemical U is ei- 
ther active or inactive: u represents the concentration of 
the chemical U in the active state and u represents that 
in the inactive state. We set the sum of u and « to 1 
(u + u = 1). An external input (S) activates U by cat- 
alyzing the process from u to u, while it is also involved 
in the synthesis of another chemical V, which acts as a 
catalyst to deactivate U. When the degradation of V is 
also included (FigfT|), the kinetics is represented as 



du 5(1 — u) 
dt r\ 



dv 
~dt 



S-v 



(2) 



We rescale time as t = t/rj in all the following discussions, 
and by setting jl = fj,/r], eq.([2]) is re- written as 



du dv 
-~ = 5(1 -u) -uv, — ~ 
dt dt 



(3) 



This equation has a stable fixed point v* — 5, u* = 
u* = 0.5, which is independent of S. When the input 
5 increases (decreases), the concentration of U in the 
active (inactive) state increases; therefore, u increases 
(decreases) from its steady state value (u* — u* = 0.5) 
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Fig. 1: (Color online) (A) schematic of a reaction system 
represented by eq.©. The area enclosed by the dashed line 
contains one element. (B) A typical example of a response 
that eq.([2| obeys; u (red x point, plotted on the vertical 
axis to the left) represents adaptive dynamics and v (blue 
solid line, plotted on the vertical axis to the right) represents 
relaxation according to the change in S (cyan dotted line, 
plotted on the vertical axis to the right). 



components in an intracellular reaction, the interaction 
represents signal transduction in an intracellular reaction 
network. To choose a specific type of interaction, we take 
the following two facts into account. First, in a biologi- 
cal reaction, input-output relations often take a sigmoidal 
form with some threshold, as in an enzyme reaction in 
which relations are modeled by a Michaelis-Menten type 
or Hill function. Further, cell-cell interactions can often 
be represented by such threshold functions, e.g., in bind- 
ing at receptors or in electric interaction in nerve cells. 
Hence, we adopt a such sigmoidal reaction to represent 
interactions between elements. The input 5 is transmit- 
ted according to the concentration of u in the active state. 
When the balanced state (it = u = u*{= 0.5)) is set to 
the threshold, the signal is transmitted depending on the 
sign of u — u* . Hence, each element receives a signal 
as a function of Y) ^ C'ij [tanh{K(it,,(i) — it*)}], where CV/ 
represents an interaction factor from element j to i, and 
the element of the interaction matrix C'ij is 0, 1, or — 1 
depending on whether the signal is non-existant, excita- 
tory, or inhibitory. We exclude self-catalysis and define 
Cu = 0. Second, catalytic reactions often occur in suc- 
cession to relay signal transduction, and therefore, the 
effects of the components are cascaded. Note also that 
the change in the signal concentration often occurs on the 
" logarithmic scale" . Without loss of generality, we can 
set 5 = 1 when there is no input from other elements, 
5 = S max > 1 if the received signal is excitatory, and 
5 = S m in = ^ ^ is inhibitory; thus, we obtain 



according to the change in 5, but after some time, it 
regains the original value it* as v relaxes to the value 
5. The dynamics represented by eq.Q have two time 
scales. The time scales are defined by the reciprocals of 
the two eigenvalues of the dynamics. The eigenvalues are 
calculated by linear stability analysis around the steady 
state. The first time scale is given by r s = 1/ (25) , which 
corresponds to the response time to the change in 5, and 
the second is given by r a = /i, which is the relaxation 
time to the original value (i.e., adaptation). Here, we 
postulate faster response and slower adaptation so that 
t s = 1/(25) < r a = /L See d, [l5[ for the relevance of 
this condition to biological adaptation. 

Each element in our coupled system follows the above 
reaction mechanism that involves two chemicals (it, v) 
and proceeds according to eq.([3]). Now, we consider in- 
teractions between such elements j and i. We assume 
that the interaction between these elements can be mod- 
eled using the input term 5, i.e., an output Uj from an 
element j triggers an input to some other element i (Si), 
and thus, the process u\, u%, •••—>• Si exists; this process 
is represented by 



Si = h(ui,u 2 , ■ ■■). 



(4) 



If each element i represents a cell, this interaction rep- 
resents cell-cell communication through signal molecules, 
whereas if each element i represents a set of chemical 



Si = exp 



Jq. C V { tanh { K K ;(*) ~ u*) }} 



(5) 



where iVj is the number of inputs that the element i re- 
ceives, i.e., Ni = Y]^—i \Cij\, and thus, comparison with 
cases involving different networks with a different num- 
ber of paths is straightforward. With this form, S max is 
given by the coupling strength a as S max — e a . Without 
loss of generality, we can set S max = 10 and a = In 10 in 
order to limit the range of 5 to 10 _1 — 10 +1 ; then, the 
form of the interaction is 



Cy[tanh{re(«3 (*)—«*)}] 



(6) 



We use k = 50 unless otherwise mentioned. The response 
time scale of each element is given by t s — 1/(25) and 
that for adaptation by t q = jl. Setting /2 = 10, we can 
use t s = 0.05 — 5 and T a = 10 in order to satisfy the con- 
dition for fast response and slow adaptation (r s < r Q ). 
There is no qualitative change in our results as long as 
the threshold-type interaction is chosen and t s < r a irre- 
spective of the parameter value. Even if the exponential 
coupling form as in eq.© is not adopted, most of the 
behaviors to be reported here are obtained. 

In the present model, each element communicates with 
other elements according to its state, i.e., depending on 
whether it is active (it > it*) or inactive [u < it*). If 
Cij > 0, an active (inactive) element j activates (in- 
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Fig. 2: (Color online) (A) interaction function represented 
by eq.([6}. The red solid line represents a positive interaction 
and the blue dotted line represents a negative interaction. 
(B) conceptual scheme of a coupled system. Each numbered 
circle represents an element with (w, v) adaptive dynamics. 
The bold red and blue dotted arrows represent positive and 
negative interactions, respectively. 



hibits) the following element i, respectively, and the con- 
verse is true for CV, < (FigHJ). We use the above nor- 
malization (summation A^) so that elements receiving a 
positive signal have identical input values and identical 
response time scales, whereas those receiving a negative 
signal have smaller input values (equal for all such sig- 
nals) corresponding to a longer time scale. 

III. RESULTS FOR TWO-ELEMENT CASE 

First, we present results for the case involving two el- 
ements (1 and 2). When one of the Cy values was 0, the 
system simply converged to the fixed point with u* = 0.5 
and v* = 1. When one of the elements does not receive 
any signal, it returns to the original fixed point, and then, 
the other element does so as well. In addition to such a 
trivial case, there are three cases for different C12 and 
C21: positive-positive, positive-negative, and negative- 
negative interactions. Depending on the type of inter- 
action, two types of behaviors, perfect adaptation and 
oscillation, are observed (Fig[3]). In perfect adaptation, 
each element eventually shows perfect adaptation, and 
therefore, the overall system also adapts perfectly and 
regains the original steady state. In the case of oscilla- 
tion, each element cannot adapt to the steady state and 




time 

Fig. 3: (Color online) Time series of u and v for a two-element 
system. A (C12 x C21 < 0) corresponds to perfect adaptation 
behavior. B (C12 > 0, C21 > 0) and C (C12 < 0, C21 < 0) 
correspond to oscillations in phase and in anti phase, respec- 
tively. In each picture, the red dashed line shows the dynam- 
ics of Mi; the blue solid line, the dynamics of 112; the magenta 
dotted line, the dynamics of vv, and the cyan chained line, 
the dynamics of V2- 



repeats the cycle of response and adaptation. The out- 
put of element 1 causes a change in the input of element 
2 according to the interaction function ([5]), and element 
2 shows response and adaptation. This output of ele- 
ment 2 in turn causes a change in Si. Thus, response 
and adaptation are repeated continuously. 

If one interaction coefficient is positive and the other is 
negative (C12C21 < 0), the system exhibits perfect adap- 
tation. Without loss of generality, we choose C21 > 
and C12 < 0. When element 1 is activated (ui > u*), ele- 
ment 2 is also activated by positive interaction (C'21 > 0). 
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On the other hand, negative interaction with element 2 
(C12 < 0) inactivates element 1. Since the inputs received 
by the elements have opposite signs, excitatory output 
is not sustained. The system then undergoes damped 
oscillations and eventually exhibits perfect adaptation. 
In contrast, the system shows a sustained oscillation if 
C12C21 > 0. Depending on the sign of C12, elements 1 
and 2 oscillate in phase (if both C12 and C21 are posi- 
tive) or in anti phase (if both C12 and C21 are negative) 
between the active and inactive states. The period of the 
limit cycle is ~ r a . 

To understand the different behaviors for different net- 
works, we performed linear stability analysis around the 
steady state (m — u* = 0.5 and Vi = v* = 1). Setting 
m = u* + Sui and Vi = v* + 5vi, the linearlized form of 
eq.© is given by 



ddui 
dt 



-28uj 



1 



-Svi 



Ka 
~~2 



dSvi 1 

— = OVi 

dt [i 



2fl^ J J 



(7) 



Defining a = (-^^C^C^i, the four eigenvalues A in 
the two-element case are given as solutions of the char- 
acteristic equation 



(A + 4) 2 (A + 2) 2 



a\ 2 = 0. 



(8) 



When C12C21 < 0, as expected, the real parts of all four 
eigenvalues are negative, and thus, the original steady 
state is stable and perfect adaptation is guaranteed. 
In contrast to a single-element case, these eigenvalues 
are complex, which leads to damped oscillation. When 
C12C21 > 0, two real eigenvalues are positive and the 
other two are negative, implying that the original steady 
state is an unstable focus, leading to a limit cycle attrac- 
tor. 

IV. RESULTS FOR THREE-ELEMENT CASE 

For a three-element system, there are six directed in- 
teractions between elements that are given by Cy , which 
takes the value 1,-1, or 0. Hence, the total number of 
possible networks is 3 6 = 729. By disregarding the cases 
that are identical due to symmetry and those in which 
one or two elements are disconnected, we obtain 78 possi- 
ble networks. For example, if the third element is driven 
by only one of the other two elements and does not trans- 
mit to other elements, this obviously results in the same 
behavior as the two-element case, and hence, such net- 
works are excluded. We have studied all of these 78 cases 
to find three distinct behaviors: (i) perfect adaptation, 

(ii) simple oscillation (as in the two-element case), and 

(iii) coexistence of rapid burst-type oscillation and a slow 
limit cycle (Fig[4]). Among the 78 systems, 7 networks 
show perfect adaptation and 46 show simple oscillation. 
These two cases are basically understood as straightfor- 
ward extensions of the corresponding two-element cases. 
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Fig. 4: (Color online) Time series of u (left) and v (right) 
for the three-element system. A (C12 = — 1, C13 = — 1, C21 = 
1, C23 = —1, C31 = 1, C32 = 1) corresponds to perfect adapta- 
tion behavior; B (C12 = —1, C13 = 0, C21 = 0, C23 = 1, C31 = 
— 1,C32 = 1) corresponds to an example of simple oscillation. 

C (C12 = — l,Cl3 = — 1, C2I = — 1, C23 = l,C3l = 1,C32 = 

0) and D (C12 = -l,Ci3 = -l,C 2 i = 0, C23 = 1, C31 = 
1, C32 = —1) correspond to a slow cycle and E (same network 
as in D) corresponds to a rapid oscillation. Note the differ- 
ence in the scales of the abscissa for different figures. In each 
picture, the red dotted line shows the dynamics of element 1; 
the cyan solid line, the dynamics of element 2; and the blue 
dashed line, the dynamics of element 3. 



6 



The remaining 25 networks show a nontrivial behavior 
that first appears in the three-element case (Fig[5]). 

There are two distinct types of dynamics depending 
on the initial conditions. In the rapid burst-type oscilla- 
tion, each Vi tends toward a fixed value v[ %x , with tiny- 
amplitude oscillations around the value, while u compo- 
nents show a large-amplitude oscillation with a period 
on the order of r s . Note that the values v( tx are distinct 
from the original fixed point value v* — 1. In the slow 
cycle, u components approach the values of the original 
fixed point, remain close to these values for some time, 
and then depart from it. This process is repeated periodi- 
cally on a time scale much longer than r a (the adaptation 
time scale). 

Now, we analyze the rapid burst-type oscillation. In 
this case, the ViS stay almost constant and no adaptation 
occurs in u^s, which oscillate on the order of the fast re- 
sponse time scale. The existence of the rapid burst-type 
oscillation is associated with the time-scale difference. As 
[1 tends to infinity, or in other words, as the time scale 
ratio T s /r a between u and v tends to zero, the influence 
of u on v through S is averaged out because u oscillates 
too fast, and hence, the oscillation in v disappears. In 
fact, we have numerically confirmed that the amplitude 
of oscillation in the w's vanishes as jl increases. 

Now, we consider this limiting behavior at T s /r a ~ 
l/jl — » 0. From the condition dv/dt = in eq. (J3J), vj tx 
is proved to be a long-time average of Si, i.e., 

v/ ix = lim - / Si(t)dt=<Si >« (9) 

(FigEJ) . Because dv/dt = 0, no adaptation occurs and 
the system is driven only by u variables with the fast 
response time scale r s ' = l/(S + v* ix ), which is obtained 
from 

^ = S(l - u) - uv flx = -(S + v flx )u + S. (10) 
dt 

Of course, the above solution should be valid only in 
the limit r s /r a — > 0, but as shown in FigO the solution 
agrees well with the present case in which T s /r a ~ (0.5 — 
0.005). In this rapid burst-type oscillation, the system 
enters a state in which the average input from all the 
other elements is cancelled; consequently, v constantly 
remains in a state that is away from the original perfect 
adaptation state. 

The rapid burst-type oscillation occurs when successive 
inputs throughout the loop of three elements change sign, 
or, in other words, there is "frustration" in the loop in- 
teractions (frustration is a term in spin glass theory|l8[). 
As a simple example, consider the case in which C%\ = 1, 
C32 = 1, and C13 = —1. When element 1 is active, 
element 2 is activated through positive interaction with 
C21 = 1 and element 3 is also activated through posi- 
tive interaction with C32 = 1; however, negative inter- 
action with C13 = —1 inactivates element 1. Hence, 
the input inconsistent with its state exists. This argu- 
ment is true even if we assume that element 1 is inactive. 
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Fig. 5: (Color online) Trajectories in the phase space of 
ui,W2,«3 (up) and of vi,V2,V3 (down) for three-element sys- 
tems of type (i) — (iii). Each orbit is drawn after neglecting 
a sufficiently long transient time. The fixed point (black x 
point) represents perfect adaptation and the limit cycle shown 
by the blue dotted line is an example of simple oscillation. An 
example of rapid oscillation is plotted as a red bold broken 
line for ui,U2,U3 and as a + point for vi,V2,V3, as the oscil- 
lation amplitude in v is negligibly small. Note that this fixed 
value of the v's differs greatly from the original fixed point 
that corresponds to perfect adaptation. The cyan solid line 
for plots of both u and v denotes a slow cycle. The results for 
(a), (b), (d), and (e) in Fig(4]are plotted. 



For any value of Ui, such opposite inputs remain, and 
hence, "frustration" among elements exists, as studied 
in spin statistical mechanics, e.g., spinglass theory. Un- 
like the case of perfect adaptation, this inconsistency in 
two-to-one anti-phase relationship cannot be eliminated, 
and hence, some outputs remain and destabilize the orig- 
inal adapted state. The v variables that cannot change 
as fast as this opposite input stay almost constant at 
v fix =< g > oc , whereas only u variables show fast oscil- 
lations. 

Next, we study the slow oscillation in which both u 
and v components show a periodic change with a time 
scale much longer than r a (the original adaptation time 
scale). As v changes slowly, the adiabatic approximation 
is valid. In particular, this approximation holds under 
the limit r s /r a ~ — > 0. Slow v variables gradually 
relax to S values with r Q , while fast u variables relax to 
the equilibrium values, defined by du/dt = 0, according 
to the values of the v variables at each instant within t s 
(FigEJ. Thus, Mi's approach the original adapted value 
and remain close to this value over the time scale r Q . For 
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Fig. 6: (Color online) Comparison of the solution v(t) (points) 
with an approximation solution of < S >oo obtained from 
dv/dt = (line). The parameters are the same as those in 
Fig(4je). The results for elements 1 (dotted line, red x), 2 
(solid line, cyan +), and 3 (dashed line, blue *) are shown. 
Notice that each approximation (line) gives a good solution 
for the behavior (points). 
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Fig. 8: (Color online) Dependence of T 3ys in the slow cy- 
cle on the strength of the interaction function (k in eq.([6j|). 
The symbols x and + represent the results corresponding to 
networks in Fig(4jc) and Fig(4jd), respectively. As the inter- 
action becomes stronger (k becomes larger), the characteristic 
time scale gets longer. The time scale of the response (r s ) is 
0.05 — 5, and the time scale of adaptation (r a ) is equal to 10. 




0.2 1 

1000 1200 1400 1600 1800 2000 
time 



Fig. 7: (Color online) Comparison of the solution u(t) (points) 
with the following solution obtained from an adiabatic ap- 
proximation: u = S/(vfi X + S) from du/dt = (line). The 
parameters are the same as those in Fig(4jd). The results for 
elements 1 (dotted line, red x), 2 (solid line, cyan +), and 3 
(dashed line, blue *) are shown. Notice that each approxima- 
tion (line) gives a good solution for the behavior (points). 



a majority of the time, ttj's are adapted to the input S 
from other elements and show only intermittent, periodic 
responses, while Uj's change gradually. By adiabatically 
eliminating the variable Ui, the motion of slow v variables 
is obtained. In FigO we have plotted the present orbit 
and the solution obtained by this adiabatic elimination; 
we find that the agreement is rather good. Note that in 
the present case, the input S varies between 0.1 and 10, 
and indeed, when S ~ 0.1, the effective response time 
t s = 1/(25") ~ T a /2. Despite this change in the time 
scale, the adiabatic approximation agrees well with our 
simulation results. 



In the slow cycle, the variables u and v undergo a repet- 
itive process in which they approach, remain nearby, and 
depart from the original fixed point (u = u* — 0.5, v = 
v* = 1); the duration of the slow cycle is much longer 
than r a . This long time scale is a result of interaction. To 
verify this claim, we studied the dependence of the period 
on several parameters in the interaction function ([5]) . We 
found that the period depends on the strength of the in- 
teraction (kct/2 in eq.© or (fTTj) ). Since r s changes with 
the value of a, we study the dependency of the period 
on the value of k. Up to some value of re, the interaction 
function h does not have a sufficiently steep threshold, 
and hence, the original fixed point is stable and the slow 
cycle state does not exist. Beyond some critical value of 
k, the slow cycle appears as a result of Hopf bifurcation. 
With a further increase in k or the interaction strength 
kct/2, the period increases even though the adaptation 
time r a remains fixed (Figj8j) . The period of burst-type 
oscillation, in contrast, remains at r s despite the change 
in k. 

For most networks that we have numerically studied, 
rapid burst- type oscillation and the slow cycle coexist. 
There exist only two attractors, depending on the ini- 
tial conditions (Fig[5]). How the two solutions are sepa- 
rated in phase space is clear when we consider the limit 
T s /T a ~ — ► 0. Under such a limit, the slow cycle 
is obtained by taking an adiabatic approximation given 
by du/dt = 0, according to the value of v at each in- 
stant. On the other hand, rapid burst- type oscillation 
is obtained in the condition dv/dt = 0. The two so- 
lutions coexist without canceling each other, since the 
rapid burst-type oscillation traces around the manifold 
spanned by Ui(i — 1,2,3) whereas the slow cycle traces 
around the manifold spanned by Vi\ the two oscillations 
occur in spaces orthogonal to each other. It is remark- 
able that two attractors whose periods differ by more 
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Fig. 9: (Color online) Basin structure of rapid oscillation (red 
x) and slow cycle (blue +) for different initial perturbations. 
The parameters are the same as those in Fig(4jd) and (e). 
The perturbation is introduced only in each of the indicated 
variables, and the remaining variables are set to the steady 
state value (u* = 0.5, v* = 1). In order to classify the attrac- 
tors, T ays is computed from the return time after a sufficiently 
long period that includes transients has been disregarded. 



than two digits coexist in the same network. 

Now, we classify the networks into classes of networks 
in which one of the following occur: (i) perfect adapta- 
tion, (ii) simple oscillation, and (iii) coexistence of rapid 
burst-type oscillation and a slow cycle; the classification 
is based on the network structure. For cases (ii) and 
(iii), we have computed the period of the system (T sys ). 
The time periods are plotted in Fig[10] as a function of a 
network index (abscissa). Since perfect adaptation corre- 
sponds to a fixed point solution, networks of type (i) are 
not included. For the parameter values chosen here, the 
time scales for a single element are t s — 1/(2S) = 0.05 — 5 
for fast response and r a = p, = 10 for slow adaptation. 
The networks with indices less than 25 belong to cate- 
gory (iii), where two types of limit cycles coexist. The 




10 20 30 40 50 60 70 
network number 

Fig. 10: (Color online) Periods T sys of three-element sys- 
tems (ordinate) plotted against the network number (ab- 
scissa). The network numbers are determined arbitrarily but 
are sorted so that No.l — 25 correspond to case (iii), where 
rapid oscillation and a slow cycle coexist, and No. 26 — 71 cor- 
respond to the simple oscillation. The networks that show 
perfect adaptation are not included here. Among the net- 
works that shows simple oscillation, No. 26 — 31 are networks 
corresponding to the shaded area in Fig llll The results for 
five different initial conditions for each network are overlaid. 
The values r s = 0.05 — 5 and r a = 10 are chosen here. The 
periods are computed from the return time after a sufficiently 
long period that includes transients has been disregarded. 



rapid oscillation has a period 0.05 $J T sys < 5, which im- 
plies that the period is within the range of r s . The slow 
cycle has a period much longer than r a : T sys > 10t o for 
the present parameter values. The coexistence of the two 
solutions on distinctly separated time scales is clearly ob- 
servable in the figure. The simple oscillation has a period 
that is slightly larger than, but on the same order as, the 
adaptation time scale r a , i.e., r a < T sys < 10r Q . For these 
networks, the simple oscillation is the unique attractor; 
as seen in the figure, the period of oscillation always lies 
between the time periods of rapid and slow oscillations. 
The classification of the type of oscillations on the basis 
of the order of T sys , r s , and r a is not influenced by the 
time scale for each element. 

To classify the cases (i),(ii), and (iii), we again per- 
formed a linear stability analysis of eq.([n]) around the 
original, perfectly adapted fixed point solution. The six 
eigenvalues A are solutions of the following characteristic 
equation: 

(A + 4) 3 (A + 2) 3 - 7 A 2 (A + I)(A + 2) - f3\ 3 = 0, (11) 



where 



P = C^f(Ci 2 C 23 C 31 + C13C32C21), (12) 
7 = (^) 2 (Ci 2 C 21 + C 23 C 32 + C 31 C 13 ). (13) 
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Note that the two indices (3 and 7 characterize a dom- 
inant network structure; (3 indicates whether the system 
has a positive or negative loop over the three elements. 
When (3 > 0, each element receives an input traversed the 
entire loop without changing its sign. Inputs to support 
and suppress activity are received through the loop when 
an element is active and inactive, respectively. Hence, 
three elements can behave in unison or in two-to-one anti- 
phase synchronization. If f3 = 0, the input through the 
loop is cancelled, while for (3 < 0, as previously men- 
tioned, the frustration already exists, and through the 
loop, each element receives an input inconsistent with 
its state. The index 7, on the other hand, represents 
the average cooperativity between the two neighboring 
elements. When 7 > 0, two elements can behave in a 
synchronized manner, while when 7 < 0, they cannot do 
so because of frustration between the two elements. 

We study the sign of the real part of the eigenvalues 
from eq. lfTTj) that determine the stability of the steady 
state and examine the dependence of the sign on (3 and 
7. Firstly, if (3 — and 7^0, there are six negative 
eigenvalues, where the perfect adaptation state is stable. 
If (3 > or (3 — and 7 > 0, there are two positive real 
eigenvalues and two pairs of complex conjugate eigenval- 
ues with negative real parts. If (3 < 0, there are four 
eigenvalues with positive real parts. The lattice region in 
Fig llll (where either both (3 and 7 are negative, or even 
if 7 is positive, the magnitude of (3 < is much larger) 
corresponds to two pairs of complex conjugate eigenval- 
ues with positive real parts. The shaded region in Fig llll 
corresponds to two real positive eigenvalues and a pair of 
complex conjugate eigenvalues with positive real parts. 

Interestingly, this classification on the basis of the 
eigenvalues of the stability matrix corresponds to the 
classification on the basis of the three behaviors (i) — 
(iii). We classify the networks that show a particular be- 
havior on the basis of the indices (3 and 7 and plot each 
type of behavior in a phase diagram in Fig llll As shown 
by the linear stability analysis, perfect adaptation occurs 
when (3 = and 7^0. Simple oscillation appears un- 
der the condition (3 > or (3 = and 7 > and also 
in the shaded region corresponding to (3 < and 7 > 0. 
The coexistence of rapid burst-type oscillation and the 
slow cycle occurs at (3 < 0, except under the conditions 
corresponding to the shaded area. 

In other words, phase (ii) appears when the system 
has no frustration (/3 > 0). Because the elements can 
cooperate even when the original fixed point is unsta- 
ble, the system belongs to the same category as that for 
a > in the two-element case, and thus, a limit cycle 
is generated. Phase (iii) appears when the system has 
frustration. There are two pairs of complex conjugate 
eigenvalues. One pair has a large positive real part and 
the corresponding eigenvector is orthogonal to the adia- 
batic space. The real part of the other pair is positive but 
has a much smaller value; the eigenvector of this pair is 
not orthogonal to the adiabatic space, and hence, in the 
adiabatic case, the instability of the fixed point solution 
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Fig. 11: (Color online) Phase diagram of the three behaviors 
in a three-element system characterized by two indices (3 (or- 
dinate axis) and 7 (abscissa axis). /3 indicates the consistency 
around the three elements and 7 shows the average coopera- 
tivity between the neighboring two elements. The definitions 
of (3 and 7 are given in the text. The green box represents 
a network that exhibits perfect adaptation behavior (i); the 
blue circles represent simple oscillation (ii); and the red xs 
represent the coexistence of a slow cycle and a rapid oscilla- 
tion (iii). 



is weak, suggesting the emergence of a slow oscillatory 
mode. Therefore, the coexistence of the two types of 
limit cycle attractors, one with a short period and the 
other with a much longer period, is a natural outcome of 
the appearance of phase (iii). 

In addition to a pair of complex conjugate eigenvalues 
with a positive real part, the shaded area in FigQj] cor- 
responds to two more positive eigenvalues that are real. 
Unlike the region for phase (iii), there is no additional 
complex conjugate pair. In fact, only a single limit cycle 
arising from the single pair of complex eigenvalues with 
a real part exists, and consequently, a simple oscillation 
is developed, as in phase (ii). The modes with two posi- 
tive real eigenvalues do not generate a different attractor. 
Even though three-body frustration exists, as indicated 
by (3 < 0, the cooperativity of two-body interactions can- 
cels the frustration. However, there still seems to be a 
difference between the behavior when [3 < and the sim- 
ple oscillation when (3 > 0. The period in this region is 
shown in Fig[l0]as the period of network (ii)'. 

V. CASES INVOLVING LARGER NUMBER 
OF ELEMENTS 

Now, we study a case involving a larger number of el- 
ements (N ^ 4). In such cases, we once again observed 
three types of behaviors; (i) perfect adaptation, (ii) sim- 
ple oscillation, and (iii) coexistence of a rapid burst-type 
oscillation and a slow cycle (Fig [T2|) . Thus far, we have 
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Fig. 12: (Color online) Examples of time series of u and v 
for four-element systems. A (C12 = 0, C13 = — l,Cu = 

0,6*21 = 0, C23 = — 1, C*24 = —1,(731 = 1,6*32 = 0, C34 = 

0,C*4i = —1,(742 = 1,C*43 = 0) corresponds to perfect adap- 
tation behavior and B (C12 = l,Ci3 = — l,(7i4 = l,C*2i = 

0, C*23 = —1,(724 = l,C3i = —1,C32 = 0, C34 = —1,(741 = 

0, C42 = 0, C43 = —1) corresponds to oscillatory behavior. C 
(Cl2 = — 1, C13 = 0,C*14 = —l,(72i = — 1, C23 = 1,(724 = 

0, C31 = 0, C32 = —1, C34 = 0, C41 = 1, C42 = — 1, C43 = — 1) 
and D (C12 = 0, C13 = 1, C*i 4 = -1, C21 = 0, C 23 = 1, C* 24 = 

0,(731 = 0, C32 = —1,(734 = —1,(741 = 1,(742 = 1, C43 = 0) 

correspond to slow-cycle behavior and E (same network as 
in (d)) corresponds to rapid oscillation. Note carefully the 
time range in each graph. In each picture, the red dotted 
line shows the dynamics of element 1; the cyan solid line, the 
dynamics of element 2; the blue dashed line, the dynamics of 
element 3; and the light-green chained line, the dynamics of 
element 4. 



not observed behaviors that deviate from the behaviors 
in these cases. 

To classify the network behaviors, we again performed 
linear stability analysis of eq.([3]) around the original, 
perfectly adapted fixed point solution. In this case, 
the dominant structures in the interaction network are 
again characterized by the coefficients of the characteris- 
tic equation, which indicate the n-loop frustration in the 
network. The coefficients are given by loop structures 
over n elements (n — 2, 3, • • • , N) containing a product 
of smaller loops as follows: 



<ij> 



£( 3 ) - (~2") 3 E C i3 C 3kCki, 
<ijk> 

£(4) = (~2~) 4 ( E CijCjkCkiCu 

<ijkl> 



(14) 
(15) 



- J2 CijCji X C kl C lk ), (16) 

<ij><kl> 

where the summation by < ij > runs over all possible 
pairs with i ^ j and that by < ijk > runs over all triplets 
with i =^ j , j ^ k and k =fi i; the pairs and triplets are 
chosen from among the N elements. 

For the case in which more than four elements are 
involved, the n-loop frustration is characterized by the 
C(n) that is determined by a product of C(rrij) (rrij ^ 

2, Ylj=i m j — n )' as follows: 
£( n ) = Ci^Ci^ ■ ■ ■ Ci n i 1 

<ii— i„> 

+ E---E 1 ] ? C( mi )...C(mu). 



mi rriM 



(_l)n+l 



(17) 



For N — 4, we again draw a phase diagram for the 
three types of behaviors, which are plotted with respect 
to the three C indices, by classifying the networks on the 
basis of the period of the attractor(s). The classification 
of the networks on the basis of the sign of L as in the 
case where N = 3 seems to be valid for the case where 
N = 4 as well (FigHSJ). 

VI. DISCUSSION AND CONCLUSION 

In this study, we have introduced a system of coupled 
adaptive elements, studied the behavior of the dynamics 
of these elements, and classified the behavior on the basis 
of its relationship with the network structure. For a sys- 
tem with two elements, two types of behaviors, perfect 
adaptation and simple oscillation, were observed. For 
a system with three elements, novel types of dynam- 
ics, namely, rapid burst-type oscillation and a slow cycle, 
were discovered in addition to the above two behaviors. 
The rapid oscillation and the slow cycle may coexist de- 
pending on the initial conditions. These dynamics are a 



11 




0.5 



CO 

i-l 



-0.5 



++ : 

+ + -U- i 
+-F : 


h + 

L++ 

f i + 


Ei-Ei— i m nm 

x .fxj 

X + |+ + ! 

* + ( 

X ! 

+ 


< + „,+ 



-0.4 -0.2 





L4 



0.2 0.4 



Fig. 13: (Color online) Classification of four-element systems 
on the basis of £(2) (abscissa; upper figure) and £(3) (ordinate; 
upper figure) as well as £(4) (abscissa; lower figure) and 
£(3) (ordinate; lower figure). The indices are defined in the 
text. The green box represents a network that exhibits per- 
fect adaptation behavior; the blue +S represent simple oscilla- 
tion; the red xs represent the coexistence of a slow cycle and 
a rapid oscillation. 1,000 networks are generated randomly 
and each network is labeled as a particular type on the basis 
of the period T sys . 



result of the existence of two distinct time scales for the 
behavior of the adaptive element, i.e., fast response and 
slow adaptation processes; hence, the dynamics are char- 
acteristic of coupled adaptive systems. The dynamics of 
each network are classified on the basis of the network 
structure (in specific, a combination of positive or neg- 
ative feedback among elements). Cooperativity among 
elements leads to simple oscillation, whereas frustration 
among elements leads to the coexistence of a rapid burst- 
ing oscillation and a slow cycle. The index for frustra- 
tion is defined by the product of signs of the interaction 
matrix elements over the entire loop. This index also 
characterizes the eigenvalues in the linear stability ma- 
trix around the steady state. 

The period of the rapid burst-type oscillation is on the 
order of the response time in a single element, implying 
the total absence of adaptation; in contrast, the period of 
the slow cycle is much larger than the original time scale 



of adaptation, where frustration among elements caused 
by interaction was important. It is remarkable that the 
attractors with such large differences in their periods gen- 
erally coexist. Frustration in a coupled dynamical sys- 
tem is also discussed in [19( , where frustration leads to a 
chaotic behavior. 

Coupled oscillators have attracted much attention and 
many general concepts as synchronization have been 
developed [l^. In oscillator networks, frustration causes 
synchronization 20] or leads to an ordered state charac- 

Synchronization in cou- 

1 m. 



terized by quasientrainment 21 1 . 
pled chaotic oscillators have also been studied [25 
However, a coupled system of adaptive elements has not 
been studied thus far. Our study here demonstrates that 
novel cooperative dynamics generally arise as a result of 
frustration in coupled adaptive systems. 

Differences in the response time scales are commonly 
observed in biological systems. The frequency of the cal- 
cium oscillations in a cell often depends on the cellular 
state, and it is proposed that this frequency is related 
to cell type differenitaion(2~ij . Neuroscience is another 
example, where two-state fluctuations in neural activ- 
ity have recently attracted much attention^, [2(| |27] : 
where a neuron spontaneously switches between the up 
state with rapid-burst firing and down state with rare fir- 
ing. Some models have been proposed to explain the me- 
chanics of two-state fluctuations by taking into account 
the recursive excitable interaction [28|, [2i|; however, the 
mechanism is not yet fully understood. Although our 
model has not been developed specifically for the neural 
system, it shows excitatory response that is common in 
the neural system. Our assertion in this context is sim- 
ple: frustration in interactions among adaptive elements 
generates bistability between burst firing and rare firing. 

The extension of the three-element system to a larger 
number of elements is straightforward. Once again, we 
observed the coexistence of burst-type oscillations and 
a slow cycle. The behaviors are classified on the basis 
of the frustration in a loop structure over n elements 
(n = 2, 3, • • ■ , N) in an interaction network; the frustra- 
tion in the interaction network is characterized by the 
coefficients in a characteristic equation in the linear sta- 
bility analysis around the steady state. These frustration 
indices are again relevant to the classification. 

The generation of long-term dynamics and differences 
in time scales is essential in determining cellular and neu- 
ral behaviors. Studies focusing on the time-scale differ- 
ence, a coupled chaotic system with a variety of time 
scales or coupled pha se-oscillators with diversification 
of the time scales [3 ll |32| have been carried out. Depend- 
ing on the network structure, such behavior is inherent to 
a coupled adaptive system in a network, despite the sim- 
plicity of the model that is utilized in this case. Here, we 
have only studied the case in which the coupling between 
identical elements is the same. In biology, elements are 
often heterogeneous and the number of elements can be 
very large. Such complex cases are currently being in- 
vestigated, and these investigations will be reported else- 
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